! VSemiletov, 08.2011
! Suggest that contour is plane and convex !
subroutine AfterReadContourData()
use Data_cfoam, only: rank, nump
use Acoustics
implicit none
include 'mpif.h'

    integer :: iPoint, i,j, ix,iy
    real(8) :: beta, r, RadVec(3), r1, L
    real(8) :: epsilon = 1.d-8

!! center of contour defintion
!    Center(1:3) = 0.0
!    do iPoint = 1,nPoints-1
!        Center(1:3) = Center(1:3) + Contour3D(1:3,iPoint)
!    end do
!    Center(1:3) = Center(1:3) / (nPoints-1.0d0)
    
! center of contour defintion
    Center(1:3) = 0.0
    L = 0d0
    do iPoint = 1,nPoints-1
        x(1:3) = Contour3D(1:3,iPoint) - Contour3D(1:3,iPoint+1)
        r = sqrt(sum(x*x))
        L = L + r
        Center(1:3) = Center(1:3) + 0.5 * (Contour3D(1:3,iPoint) + Contour3D(1:3,iPoint+1)) * r
    end do
    Center(1:3) = Center(1:3) / L ! (nPoints-1.0d0)

    if(rank.eq.0)then
        write(*,*) ' '
        write(*,'(3f25.10)') Center(1:3)
        write(*,*) ' '
    end if

! 2D basis definition

    r = 1.0d10
    ix = 0
    iy = 0
    do i = 1,nPoints-2
        do j = i+1,nPoints-1
            x(1:3) = Contour3D(1:3,i) - Center(1:3)
            y(1:3) = Contour3D(1:3,j) - Center(1:3)
            r1 = abs(sum(x(1:3)*y(1:3)))/sqrt(sum(x(1:3)*x(1:3))*sum(y(1:3)*y(1:3)))
            if(r.gt.r1)then
                ix = i
                iy = j
                r = r1
            end if
        end do
    end do

    x(1:3) = Contour3D(1:3,ix) - Center(1:3)
    y(1:3) = Contour3D(1:3,iy) - Center(1:3)

!    x(1:3) = Contour3D(1:3,1) - Center(1:3)
!    y(1:3) = Contour3D(1:3,2) - Center(1:3)
    
    r = sqrt(sum(x(1:3)*x(1:3)))
    if(r.lt.epsilon)then
        stop 'Fatal error: Contour is incorrect!'
    end if
    x(1:3) = x(1:3) / r
    beta = - sum(x(1:3)*y(1:3))

    y(1:3) = y(1:3) + beta * x(1:3)
    r = sqrt(sum(y(1:3)*y(1:3)))
    if(r.lt.epsilon)then
        stop 'Fatal error: Contour is incorrect!'
    end if
    y(1:3) = y(1:3) / r

! 2D contour definitioin
    do iPoint = 1,nPoints
        RadVec(1:3) = Contour3D(1:3,iPoint) - Center(1:3)
        Contour2D(1,iPoint) = sum(RadVec(1:3)*x(1:3))
        Contour2D(2,iPoint) = sum(RadVec(1:3)*y(1:3))
    end do
    
! Center direction
    CenterDir(1:3) = CenterDir(1:3) / sqrt(sum(CenterDir(1:3)*CenterDir(1:3)))
    
!    Center(1:3) = Center(1:3) + 0.5d0 * Length * CenterDir(1:3)
    
end subroutine AfterReadContourData


!   	type TAcousticSurf ! 
!        logical :: active=.true.  ! on/off
!        integer :: nPoints ! number of countour points
!		real(8) :: Contour3D(3,1000) ! 3D coord of convex contour
!        real(8) :: Center(3)         ! coord Center of convex contour // left disc's center
!        real(8) :: x(3),y(3)   ! 2D basis in plane of the convex contour
!        real(8) :: Contour2D(2,1000) ! 2D coord of the convex contour in 2D basis
!! input here! NOTE, Semiletov
!        ! radius-vector from th left disc's center to the right disc's center
!        real(8) :: CenterDir(3)=(/0.,0.,1.0/)
!        ! factor of contour decreasing from the right side
!        real(8) :: factor=1.
!    end type TAcousticSurf
